Empirical parameterisation and dynamical analysis of the allometric Rosenzweig-MacArthur equations

Allometric settings of population dynamics models are appealing due to their parsimonious nature and broad utility when studying system level effects. Here, we parameterise the size-scaled Rosenzweig-MacArthur differential equations to eliminate prey-mass dependency, facilitating an in depth analytic study of the equations which incorporates scaling parameters’ contributions to coexistence. We define the functional response term to match empirical findings, and examine situations where metabolic theory derivations and observation diverge. The dynamical properties of the Rosenzweig-MacArthur system, encompassing the distribution of size-abundance equilibria, the scaling of period and amplitude of population cycling, and relationships between predator and prey abundances, are consistent with empirical observation. Our parameterisation is an accurate minimal model across 15+ orders of mass magnitude.


Introduction
Allometric scaling relationships have been the subject of scrutiny and debate since the connection between organism size and its metabolic rate was first defined by Rubner in 1883 [1][2][3][4]. These models, which link some characteristic y to the size x of an organism via the power law y = ax b (where a, b are scalar constants), are appealing due to their capacity to capture a multitude of relationships despite their simplicity. Scaling laws have been used to express a variety of biological rate measures, such as metabolism, consumption, and birth or death rates [5][6][7][8].
Allometry is also utilised in modelling behavioural traits and bioenergetic characteristics, such as movement behaviour or locomotory costs [7,[9][10][11]. At broad scales, such laws have been applied to ecosystem-level properties, including predictions of organism population density and carrying capacity [12][13][14]. However, despite scaling laws' wide utility and intensive study, there has been a limited exploration of the properties of minimally constructed, size-generalised predator-prey models [15][16][17]. Many authors have examined the empirical relationship between organism and population sizes [12,13,[18][19][20]. Reported exponents fall between −1 and −1/4 depending on factors such as taxonomy or environment. The classical −3/4 value describing global size-density relationships is the direct inverse of Kleiber's 3/4 law for metabolic scaling [2,13], leading to the 'energetic equivalence' hypothesis: that is, the net energy contained within each size class is invariant [18,20]. This conjecture has been widely debated, particularly with respect to whether this invariance is cause or effect of other bioenergetic drivers [21]. However, despite disagreement over underlying mechanisms, there is broad consensus that the consistency of size-density scaling within empirical data likely reflects fundamental physical constraints [5,21]. To examine what drives limitations in macro-scaling behaviour, it is possible to use dynamical size-based models incorporating organism traits that scale across the size range [6,15,16]. This approach facilitates the investigation of critical breaks in ecosystem-level scaling laws within a global framework, and the exploration of potential impacts from changes that may affect many organisms in a similar way-for example, warming temperatures or emergent hypoxia in the oceans [22,23]. However, perturbing parameters across 15+ orders of magnitude in size poses challenges. For example, coexistence regions of size-generalised predatorprey models are dominated by scaling exponents [16].
It is more straightforward to keep model behaviour stable, thus resolving the coexistence issue, by using the 4-parameter Lotka-Volterra model, but that setting is too simple for some applications [6,17]. Alternately, a series of models may be solved piece-wise for different size classes, yet this means that they are not truly generalised. In the most comprehensive study to date, a size-based paramaterisation of the Rosenzweig-MacArthur system places restrictions on the relationships between the exponents of each parameter [16]. However, this in turn limits the types of perturbations that may be applied or investigated. Finally, there are discrepancies in the treatment of the functional response term between the theoretical and empirical literature. The theoretical literature broadly assumes that the limit of maximal consumption ties predator production to the prey's and thus scales negatively to match prey production, however, there is empirical biological evidence for positive scaling [6].
Here, we present an alternate approach of paramaterising size-based predator-prey interactions for the classical Rosenzweig-MacArthur differential equations. Under this framework, we use multiple forms of analysis to create a robust picture of parameter sensitivity and required ranges for species coexistence in the context of real-world observations. We are able to show that, despite the number of assumptions inherent within this style of modelling, the mathematical restrictions are closely related to biological observations. Finally, we describe the conditions required to create an entirely size-invariant model, and show how empirically derived parameters generate ODE solutions that match real-world size-abundance distributions.

Parameterisation of the model
We begin with the Rosenzweig-MacArthur ODEs [24] and Holling II functional response, We have variables R for resources and C for consumers. The parameters r and δ are birth and death rates respectively. Carrying capacity is given by K, interaction rate b, handling time h and the conversion efficiency �. To investigate the system across the full size range we scale the parameters by mass. Organism size (in g) is given by S R for resources and S C for consumers. We depart from [16] by constructing the functional response term in line with parameterisations used within experimental research [6,25,26]. Hence, the strictly positive parameters are expressed as where for each parameter i, the coefficients i 0 may be standardised (Appendix A in S1 Appendix), and α i denotes the scaling exponent. Next, we define the prey-predator mass ratio as ρ, where ρ > 0. We may then relabel the parameters r, h and K in terms of the consumer, With this approach we extend the results of [16] by placing no restrictions on the exponents, allowing h to be an independent term which may be matched to empirical observations. We now also follow standard practice by setting � / ρ, that is, the conversion efficiency is proportional to the prey-predator mass ratio [16]. Next, we use a standard rescaling of (1) to reduce the number of parameters and simplify analyses. We setR ¼ 1=ðbĥÞ,C = �=ðbĥÞ, m ¼K bĥ, and define u ¼ R=R and v ¼ C=C. After scaling time byr such that t ¼rs, and defining g ¼ �=ðĥrÞ and o ¼ d=r, we arrive at the new system The parameters in (4) are also all strictly positive and scale across the size range. For completeness, we provide the explicit relationship between the old and new parameters in Table 1, and Table 1. Relationship between parameters in original and rescaled Rosenzweig-MacArthur system. Here, α h = α hR + α hC . the system of equations with substituted terms is

Definition Coefficient Exponent
The expression S R = ρS C facilitates interpretability in downstream analyses. All exponent terms may be collected within S C , which we henceforth refer to as S. This provides simplified expressions within Table 1 and (5), yet we may still examine the impacts of perturbations to any one parameter. Table 2 summarises parameter exponent ranges within empirical research. Code used to generate the following results and figures is available at https://github.com/ jcmckerral/universalallometry.

Coexistence & sensitivity
The non-trivial equilibrium of interest (coexistence) is obtained by equating the right side of (4) to zero and solving for u = u � and v = v � , yielding  (4) where extreme min-max limits are given for completeness. However, there is consensus that α r , α δ ' −1/4, also verified in a substantial recent review [27]. Similarly, despite the potential range for α b and α h , in large generalised studies typically 1/2 � α b � 1 [25,26], and α h � 1/8 [6,25,26] which significantly constrain the exponent ranges in (4). The bottom section of the table therefore provides exponent values in the rescaled system based on upper/lower bounds for the most plausible generalised empirical scaling values for (2) (as determined by cited articles), i.e. those found from studies with large quantities of data, across broad taxonomic and size ranges, and/or named outliers being excluded [6,22,[25][26][27][28][29]. For there to be non-negative values for (u � , v � ), we require that

Symbol
and For the Jacobian of the right side of (4) evaluated at the equilibrium point (u � , v � ), the condition det > 0 is also fulfilled by (7b). We may express (7b) as γ/ω > 1+ 1/μ. As all parameters are strictly positive, if (7b) is satisfied, it immediately follows that (7a) is satisfied also. The inequality determines the sign of the trace of the aforementioned Jacobian, which dictates whether the system converges to a point or to a stable limit cycle, and there is a Hopf bifurcation at equality. The dynamical characteristics of the Rosenzweig-MacArthur system have been explored in depth elsewhere (such as [24,32,33] and references within). Our focus is the interplay between biological and mathematical constraints. We now interrogate the behaviour of the system and inequalities (7) and (8) using a series of complementary analyses, and discuss the mathematical implications in the context of empirical observations. Handling time. The condition from (7a) is equivalent to When investigating coexistence across the large size ranges considered in this study (>15 orders of magnitude), the coefficients on the left hand side of this inequality will have negligible impact relative to the scaling exponents of the parameters, as in log space they may translate the intercept of the defined line but cannot change its slope. Therefore, we instead examine the exponents given in (9) and observe that if α δ + α h < 0 the condition may fail for small organisms. As α δ is tightly constrained ( Table 2, and comprehensively reviewed in [27]), to better understand the contribution of handling time to coexistence constraints, we now investigate system behaviour when varying the scaling of h, under the assumption that other exponents are either fixed or undergo only minor perturbations; their relative contributions to coexistence properties are explored in later sections. Handling time's classical null model from Yodzis & Innes [15] based on metabolic theory is , orĥ / S 1=4 as derived in [25]. A null model derivation assumes the maximal consumption rate of the predator-the inverse of handling time-scales with metabolic demand S 3=4 C , and the per-prey metabolic demand is therefore S 3=4 C S À 1 R . This matches the assumptions of [16,17], where the birth rate of the predator in the presence of unlimited resources is assumed to scale with the birth rate of prey, though it should be noted that other interpretations of the same model either do not normalise against prey mass e.g. [31] or do so implicitly e.g. [6]. However, consideration of physiological traits' nuanced impact on maximal consumption, and equivalently handling time, has since led to a departure from the traditional metabolic framework. Attacking, killing, then eating and digesting prey all impact the parameter h [26] and the prey's contribution to the process should be incorporated [25,26]. Handling time has therefore been recast to the more biologically representative form we use in this study: h / S a hR R S a hC C , where typically 0 � α hR � 1 and −1 � α hC � 0 [25,26]. The exponents α hR , α hC have been empirically determined in several reviews and display considerable variability [6,25,26,34]. We now discuss the implications this variability has for coexistence under the inequalities in (7a) and (7b). Two of the three reviews listed above conclude that the predator-prey components of handling time scale more gently (whether positive, or negative) than null models predict. However, the resultant exponent forĥ is positive ('1/3) for arthropod functional responses [25], but negative ('−1/8) when examining broader taxonomic groups [26]. Most of the organisms in [26] display negative scaling forĥ in taxa-specific breakdowns due to gentler scaling of the resource exponent. Only α hC is assessed in [6], which is calculated across a wide range of taxa for 2D and 3D environments and for a larger mass range than in [26]; for the 2D and 3D case α hC ' −1.1. The authors account for the steeper scaling relative to metabolic expectation by noting that feeding is an active process scaling with maximal rather than basal metabolism [6].
To calculate the scaling ofĥ from the empirical assessment in [6], we use the assumption α hR = 1, which is the parameter's upper limit. This implies the exponent α h � −0.1, and thatĥ scales below the value of 1/4 assumed by previous theoretical work on the model. Conceptually, this indicates that the parameter may be constrained by physical processes rather than a bioenergetic flux balance. Note that despite the phenomenological formulation of the functional response predator production is implicitly constrained by the prey density. We next examine coexistence condition (7b), which may be expressed as γ/ω > 1 + 1/μ meaning ln(γ/ω) > ln(1 + μ −1 ); if we then substitute the original parameters we arrive at the inequality where c 1 , c 2 are constants derived from the coefficients. Considering (10) together with the empirical behaviour ofĥ, if α h ' −1/5 or less across the size range, smaller organisms may violate this condition (Fig 1). Similar to the behaviour of the first inequality (9) in this section, significantly perturbing the coefficients does not qualitatively change this behaviour. For example, a change of 3 or 4  Table 2, here α K and α b have less effect than α h ; we thus set α K = −3/4 and α b = 1/2. https://doi.org/10.1371/journal.pone.0279838.g001 orders of magnitude in the coefficient ratios will result in a vertical translation of the line ln(c 1 ) of 3 to 4, which may impact results if one were to consider the point of intersection with f(S) for restricted size ranges (e.g. within restricted taxonomic groups), but do not affect the global behaviours of the complete size distribution considered in this study. In addition, even if the coefficients were to vary by an arbitrarily large amount, should the variation remain proportionally consistent, the ratios would remain constant; that is, the position of the coexistence line would not change. Conversely, Fig 1 shows that a relatively small perturbation of 0.6 in α h alters the value of the term f(S) by over 10 across the size range, with a particularly significant impact on the smallest organisms. Therefore, provided that α h ' −α δ then all sizes will fulfil the condition, allowing coexistence across the full span of the model. There is empirical support for these observations. The taxonomic group breakdowns in [26] indicate that smaller taxa may display positive scaling forĥ as concluded in [25] and the strongly negative scaling is observed in macro-organisms, particularly vertebrates. The notable exception of unicellular marine organisms (α h ' −1/3) has a very small sample size; further experimental studies and especially those which generate larger data sets may potentially reach alternate conclusions and we leave this for future work.
Next, we examine the contributions of the carrying capacity and interaction rate parameters K and b to coexistence and dynamical properties of the model.
Carrying capacity and scaling of population cycling. As the right side of (7b) is equal to u � , the coexistence condition may also be expressed as where u � / S a d þa h . If we combine (7b) and (8), we obtain Together, (11) and (12) indicate the carrying capacity must be sufficiently high for a sustainable prey population, and that predator attack rates must be high enough to compensate for mortality across all sizes. Assuming reasonable values for α K , α b such as those given in the body and caption of Table 2 respectively, these inequalities will generally hold. Whilst the lower bound of α b ' 1/2, estimates of the 'universal' value from comprehensive reviews suggest a number in the range 0.6 < α b < 0.9 [25,26]. This concurs with the properties of inequality (12), which-even if worst-case values are set for K and h-allows for coexistence across the full size range provided α b does not exceed 0.9 by a significant margin, or equivalently, α μ does not exceed '0.2. Smaller (or even negative) values of α b serve to improve coexistence properties, suggesting that generalised empirical values for α b conform well to the mathematical properties of the model, sitting at an upper bound of '0.9. Under the original system (1) and assuming coexistence, resource equilibria will scale with size as R � / S a d À a b and consumer equilibria will scale as C � / S a r À a b . Despite their importance for the coexistence domain, the carrying capacity and half saturation do not meaningfully impact equilibria abundances in allometric paramaterisations. However, they do impact some properties of the limit cycle. Allometric settings of the Rosenzweig-MacArthur system usually result in oscillating solutions due to size-scaled parameter values relative to the constraints in (8) [28]. Using data extracted from the literature for a broad range of taxa (Appendix A in S1 Appendix , Fig 2b), we find stronger empirical support than in previous work for a S 1/4 scaling signal for the period τ t [28], noting that τ t indicates that we are considering period under the timescale of variable t. Expressions for theoretical scaling whereby t t / ðS a d þa r Þ À 1=2 have previously been given in the literature [15,28] and an exact equation for τ t in [16], but to our knowledge a complete derivation for the equation in [16] has not previously been published. For the reader's interest we provide one in Appendix B in S1 Appendix, where we show that, under the assumption that S C / S R , ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi �=ðĥdÞ þ 1 matching the findings of [16]. We find the theory agrees well with observed values, as our empirical data has a relationship t À 1 t / S À 0:2 (Fig 2b). Whilst multiple factors are proposed to affect cycle behaviour, these factors are made up of internal and external effects [17,35,36]. In particular, cycle periods are proposed to be predominantly driven by internal physiology such as maternal generation time, whereas cycle amplitude is thought to be dominated by external effects such as predation or other environmental influences [37]. With some caveats due to a temperature dependency and clear contribution of life history traits, at a broad scale maternal effects (i.e. generation times) are suggested to mirror the inverse of the maximal growth rates r, that is, the period should fall close to S 1/4 as predicted by the model [35,38,39]. The empirical evidence for 1/4 power scaling is particularly strong for mammals [28,35,37,39]. We note that most prior work has been restricted to specific taxonomic groups or across smaller size ranges than we consider here, as we include prokaryotes, protists, invertebrates as well as herbivorous and carnivorous mammals in our dataset (Appendix A in S1 Appendix). The notable exception of [38], with a generously sized and wide-ranging dataset, finds a slightly steeper scaling of τ t / S 0.31 across aggregated taxa groups. However, the within-group exponents (considering insects, zooplankton/protists, and vertebrates independently) are shallower at '0.2 on average, which may suggest that there were challenges associated with aggregating data across different studies. That being said, period scaling of 0.2 to 0.3 across various taxa groups is within expected empirical variability, especially when studies have limited data, and we expect that the generation and synthesis of large, broad datasets may shed light on the validity of the τ t / S 1/4 Rosenzweig-MacArthur model prediction in future.
With respect to the cycle amplitudes, previous work has found that the ratio of maximum to minimum densities is size-invariant [28], indicating that the oscillation amplitude decreases with increasing size. Further qualitative support that this mathematical behaviour is aligned with features of real-world biological systems is given by the fact that log-transformed sizedensity relationships demonstrate near-constant variance across 15+ orders of magnitude [13,27]. Under our paramaterisation the log-scaled oscillations are relatively sinusoidal and may be constrained between one to three orders of magnitude (Fig 2a). Hence, they are more realistic than allometric Lotka-Volterra dynamics which push predator populations to unreasonably low levels with fluctuations exceeding 15 orders of magnitude [17]. Unfortunately early efforts to find analytic approximations of the oscillation amplitude of the Rosenzweig-MacArthur system have not been generalised [49]. Recent results have only been derived for specific-and restricted-parameter values [50]. However, the rescaled system (4) provides scope for us to examine the effects of perturbations in a simplified manner. In Fig 2c, we show through simulation that perturbations to all parameters impact the magnitude of the fluctuation of the limit cycle. However, unless these perturbations are applied to α r or α δ , the oscillation amplitude will, within a small noise factor, remain invariant with respect to the mean population density, which reflects empirical findings [28]. Ecological theory and observation suggest that oscillation amplitudes are impacted by a myriad of environmental factors, reflecting the fact that the fluctuation size in the system may be impacted by changes to any parameter. Thus, the qualitative behaviour of the Rosenzweig-MacArthur's limit cycle better reflects known biology than the Lotka-Volterra system [17,37]. We next use a sensitivity analysis to assess the system's robustness to different forms of perturbation.
Sensitivity. A local sensitivity analysis provides a first-order approximation of the relative impact of changing parameters on the solutions of (4) near the system's equilibria. We adhered to the methodology described in [51]. In order to check how sensitive the system, _ x, is to small changes in parameters, λ i , we construct a sensitivity function, S(t), such that and x(t, λ) is a solution of _ x. Next, we characterise the solution to the sensitivity equation given by _ SðtÞ ¼ Aðt; l 0 ÞSðtÞ þ Bðt; l 0 Þ; Sðt 0 Þ ¼ 0: Applying (14) to (4), A is the Jacobian of (4) with respect to variables u and v, and B is the Jacobian of (4) with respect to parameters μ, γ, and ω, both of which are evaluated at nominal parameter values. After setting initial conditions u 0 and v 0 , we obtain numerical solutions for (14). Fig 3 shows the trajectories for two initial conditions, noting that these are not the trajectories of the rescaled Rosenzweig-MacArthur system (4), but instead show the relative impacts of perturbing different parameters on the long term evolution of its solutions. Firstly, we show u 0 and v 0 in the neighborhood of u � , v � respectively (Fig 3a and 3b), and secondly for u 0 , v 0 an order of magnitude greater/smaller than u � , v � respectively (Fig 3c and 3d). The qualitative behaviour remains the same in both cases. For an initial state near the equilibrium (Fig 3a and  3b), there is monotonic behaviour as the system converges to the limit cycle. In the case of Fig  3c and 3d the limit cycle emerges after some critical time t c . This analysis reflects the globally stable nature of the Rosenzweig-MacArthur equations and shows that the system is least sensitive to μ, providing further support that perturbations to r and δ have the largest potential impacts on its dynamics. The qualitative behaviour is similar for other nominal parameter values, provided they are not set on the other side of the bifurcation boundary. The exponents with the least empirical variation-by a significant margin-are α r and α δ , mirroring the mathematical constraints. That is, the dynamical behaviour of the system is relatively robust to perturbing the functional response parameters displaying the highest empirical variance.

Applications
The rescaled parameter definitions in Table 1 may be used to determine a size-invariant system by balancing the exponents. This is desirable as it is straightforward to set coexistence for an arbitrary size range: the rescaled system remains fixed, and transforming back to the original parameter space merely translates and stretches (or squashes) the dynamics in the base coordinates. This facilitates analytic study of the equations as it greatly simplifies the task at hand: only one equation (the size-invariant version of (4), which has constant parameters) needs be studied to understand the dynamical behaviour of (1) for any sized predator or prey. This is certainly mathematically expedient and we note that the majority of allometric studies to date have used this approach. To determine the scaling values necessary, the definition of o ¼ d=r indicates that the size scaling of r must match δ, both of which consistently display an exponent of −1/4. It immediately follows that α h = 1/4 as g ¼ �=ĥr. Assuming that carrying capacity scales with a −3/4 exponent, and using the definition m ¼Kĥb it follows that α b = 1/2. Indeed, previous analytic research into scaling dynamics of the Rosenzweig-MacArthur system places limitations around scaling values using similar principles, which reduces the complexity of the analysis but potentially decreases the model's applicability to real-world systems [16]. In particular, given the tenuous empirical support forĥ / S 1=4 , it may be more biologically sound to assign a value where α h � 0; furthermore, our results from Sections 3.1.1 and 3.1.2 indicate that coexistence across the domain should be feasible under this condition. We find that it is still possible to generate a full size-abundance distribution when α h � 0. A recent size-density scaling analysis indicates that the relationship follows N � / S −1 rather than the canonical N � / S −3/4 (where N � is population density) [27], and we will show that we are able to reproduce this when choosing empirical scaling values across all parameters.
Having assessed coexistence properties in previous sections, we now note a possible extension to parameter b (interaction rate) and also consider whether our treatment of � (conversion efficiency) is reasonable for the model. Whilst a static value in the range 0.6 < α b < 0.9 is accepted as a reasonable generalised exponent based on major reviews [25,26], a limitation of our treatment of b is that we place no restrictions on the interactions between a predator and any arbitrary-sized prey. Predator-prey interaction processes are complex and increasing evidence suggests they follow a 'hump-shaped' curve with the predator-prey mass ratio [25]. A natural extension to our model would be to introduce a term reliant on ρ to the parameter b, and f(ρ, b 0 ) is a function assigning probability of prey capture based on the prey-predator size ratio. While a form for f(ρ) has been proposed [16], it would be possible to use a function encoding a broader range of life history traits for the tradeoff of introducing more parameters. Relevant processes to consider may include habitat effects on foraging, prey refuges, and optimal size ratios, resulting in further constraints on coexistence [25,26,34]. However, as this would introduce additional complexity, for the purpose of this study we assign a constant value for the coefficent, meaning that interaction rates are consistent across different prey-predator mass ratios and the exponent α b = 2/3.
Our final consideration is the contribution of ρ to the conversion efficiency �. The empirical distribution of ρ is approximately lognormal peaking at '0.02 [29]. Equilibria population ratios do not follow a 1:1 relationship with the size ratio of prey and predator when using the relationship � = ρ (solid white line, Fig 4). For a fixed predator size and increasing prey size, organisms become less efficient at converting biomass. However, this result does not align with observed data. A review of 15,000+ predator-prey pairs concludes that size differences between predator and prey has an upper limit, potentially due to inefficiencies when the size discrepancy becomes too extreme [52]. Furthermore, the larger the predator, the more generalist its feeding strategies [25]; increases in prey biomass-which could indicate predators feeding on smaller prey-do not translate to a proportionate increase in predator biomass [53]. We therefore apply the assumption that energetic reward (and biomass conversion) for predator effort declines as the size difference increases, and that the scaling of ρ with equilibria population ratios is superlinear. We can implicitly capture the result by assigning a function � = aρ ψ , where ψ is a scalar. For simplicity, we set a to 1, and in Fig 4, we assess the predatorprey population ratios for varying ψ.
A value of ψ > 1 increases the difference between the predator-prey populations; ψ < 1 reduces it. More sophisticated functional forms may include favourable size ratios or introduce a size dependency to the value of �. However, there is limited empirical research on scaling properties of � [13,19,30]. A theoretical investigation of optimal predator-prey size ratios together with more complex functional response formulations reflecting alternate foraging/ feeding strategies may yield interesting results. We leave this question open for future work.
To generate the size-abundance distribution shown in Fig 5, we use empirically motivated values where ψ ' 1.3, α r = α δ = −1/4, α h = −0.1, α K = −3/4 and α b = 2/3. Coefficients are standardised to boundary values (Appendix A in S1 Appendix). The model's distribution scales to −0.95, matching the value observed in [27]. The inset (generated from identical parameter values) shows the predator-prey density relationship, scaling at 0.76, close to the '3/4 findings in [53]. Our choice of ρ = 0.02 was motivated by the fact it is the most commonly observed preypredator size ratio [52]. Changing the value of ρ will slowly perturb the predator-prey density slope due to ψ influencing the maximum and minimum values in the limit cycles in a nonlinear fashion (also seen in Fig 2c).
We propose that the interplay between scaling of ρ and predator-prey density slopes will have mathematically rich behaviours, especially as there is not yet a general solution for limit cycle amplitude in the Rosenzweig-MacArthur system [49,50], but it is beyond the scope of this work to interrogate this question in detail. Our empirically-driven paramaterisation for We show examples of five predator sizes: 1E-6 (�), 1E-4(×), 1E-2 (+), 1E0 (4), and 1E2g ( � ). For each ρ = 0.02. The slope within each symbol group '0.76. That is, increasing prey density does not result in a 1:1 increase in predator density. The scaling relationship is sublinear, matching the observations of [53].
https://doi.org/10.1371/journal.pone.0279838.g005 our model generates a theoretical size-abundance distribution matching the largest data study to date [27], and we note that to generate Damuth's law [12], simply changing the scaling of one parameter b such that α b = 1/2 results in exponents of −0.81 and 3/4 for the size-abundance and predator-prey density scaling respectively. This sensitivity of the slope to interaction rate is reflected in biology. Variability in hunting and feeding strategies is a plausible reason that taxonomy-specific size-abundance distributions have exponents ranging from -1 to -1/4 as interaction rates are a critical driver for obtaining the energy for reproduction [12,13,[18][19][20]. The functional response literature implies that the scaling of b is shallower for terrestrial endotherms than across broad taxonomic groups [26], which is precisely reflected by the model's size-abundance scaling behaviour: Damuth's original research focussed on mammals [12], yet Hatton et al.'s incorporates the full spectrum of eukaryotes [27]. We therefore suggest that the Rosenzweig-MacArthur system is suitable not only for large scale studies such as this, but that our results may be modified for subsets of taxa by using parameter scaling specific to those organisms, creating biologically meaningful models that may be used for prediction.

Conclusions
Here, we investigate the links between empirical and theoretical allometric literature. The resource size dependency is eliminated from the system by explicitly encoding the prey-predator mass ratio, ρ. This simplifies analyses, allowing us to extend previous theoretical work by removing all restrictions on parameter relationships when examining properties of coexistence and macro scaling behaviours. Taxon-specific applications (which are far more sensitive to perturbations in coefficients due to the small size ranges under consideration) are likely to fall within the noise factor of the large size domains considered in this paper, however our methodology may still provide a parsimonious base for customising the equations in those settings. This may be useful for food web or trophic modelling, especially to mitigate against overfitting challenges [54]. Interrogating the model behaviour through three separate analyses provides a level of robustness to the finding that the mathematical constraints complement empirical observation. Contrary to most previous studies, we use an empirically determined parameterisation of the functional response term. Our results suggest that the standard approach of setting exponents based on metabolic theory may need to be reassessed [16]. The handling time parameter shows the greatest departure from those assumptions, and the highest variance, which is consistent with the massive trait variation in foraging strategies. Nevertheless, we find that results generated from an empirical setting agree with results in recent reviews of sizeabundance scaling. This work may be extended in several ways. Firstly, one could incorporate temperature effects, for example after [22,55], which may further stabilise the model by reducing the interaction strengths [30]. Secondly, additional empirical data on functional responses at the size extrema could more accurately define the scaling ofĥ and b. Type I, Type III, or generalised functional responses may also be examined, although we note that system behaviour usually remains qualitatively similar over large size ranges [6,31].
The broad limitation of allometry is that a generalist strategy can be a poor predictor of taxon-specific outcomes. Challenges to the framework arise not only from biological differences but also from physical or spatial processes, such as prey patchiness or heterogeneous habitat distribution [56]. Thus, care over interpretation and the applicability of results must be taken, particularly at the size limits in either direction. For example, prokaryotic reproduction rates fall between minutes and millenia [57,58]. Furthermore, large organisms such as whales play a critical role in nutrient recycling; assuming a single species may be defined as a resource or consumer alone does not account for the intrinsic complexities within natural environments [59]. However, despite these caveats, allometric approaches have been found to outperform those explicitly encoding organisms' individual and life-history traits when investigating a system's macro properties [25]. Classical population dynamics models remain a powerful tool in ecology, and the consistency across many allometric laws suggest self-organising processes we are yet to unravel. We propose that systematically assessing where theoretical and empirical properties of allometric modelling diverge may assist in identifying plausible mechanisms governing these phenomena.